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iH ' Abstract 

We report the results of a new accurate evaluation of light nuclei yields in primor- 

0\ \ dial nucleosynthesis. All radiative effects, finite nucleon mass, thermal and plasma 

corrections are included in the proton to neutron conversion rates. The relic densi- 

J^ ■ ties of 4 He, D and 7 Li have been numerically obtained via a new updated version of 

the standard BBN code. In particular the theoretical uncertainty on A He is reduced 

O ' to the order of 0.1 %. 
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1 Introduction 

Big Bang Nucleosynthesis (BBN) still represents one of the key subject of modern cos- 
mology even if its clear understanding traces back to over 25 years ago Jl|. The reason for 
this relies on the fact that BBN is one of the most powerful tools to study fundamental 
interactions, since light nuclei abundances are crucially depending on many elementary 
particle properties. As a well-known example, the 4 He abundance is strongly affected 
by the number of effective neutrino degrees of freedom, but others fascinating phenom- 
ena such as neutrino degeneracy or oscillation phenomena can be studied too, using the 
universe few seconds after the bang as a laboratory. 

In the recent years, the experimental accuracy of the light primordial nuclei abun- 
dances, mainly the one of A He, underwent a sort of revolution. The qualitative results 
of not too many years ago, suggesting that the 4 H e mass fraction I4 was of the order 
of 0.25, recently turned in measurements with accuracies of the order of one percent. A 
similar good improvement has been obtained in both Deuterium (D) and 7 Li abundances 
Y" 2 = D/H and Y 7 = 7 Li/H. In particular for D, measurements in distant Quasars 
Absorption line Systems (QAS) now represent a reliable estimate of the primordial value 
for Y 2 , which is only lowered by subsequent stellar processing. Paradoxically, the refine- 
ment of these experimental techniques, due to the uncertainties in the models describing 
stellar activity, is at the basis of large discrepancies between different set of results. Such 
discrepancies are possibly of systematic origin, or may reveal new aspects of cosmological 
evolution of the universe. The observations of I4 from regression to zero metallicity in 
blue compact galaxies in two independent surveys still produce two incompatible results, 
a low value 01, 



Y} 1) = 0.234±0.002±0.005 , (1.1) 



and a significantly higher one 



Y} h) = 0.243±0.003 . (1.2) 



A similar situation occurs in D measurements, where observations in different QAS, both 
at red shift larger than 3, give two results at bias for one order of magnitude |jj, [5] 

F 2 (0 = (3.4±0.3) 1(T 5 , (1.3) 



Y 2 {h) = (1.9±0.4) 1(T 4 . (1.4) 

For the 7 Li abundance, the almost constant Spite plateau observed in the halo of POP II 

stars H EU 

Y 7 {1) = (1.6±0.36) 1(T 10 , (1.5) 

is generally considered a reliable estimate of primordial abundance. Nevertheless, the 
observation of stars similar to the ones contributing to the Spite plateau, but with no 
traces of 7 Li || [5], seems to imply the presence of a depletion mechanism. A recent 
analysis based on a sample of 41 stars does not find any evidence of depletion mechanism 
or post-BBN creation and yields the primordial abundance || 

YJ h) = (1.73±0.21) 1(T 10 . (1.6) 

A brief summary of the complete experimental situation on primordial abundances can 
be found in Ref. [pT0f| . 

Probably, future measurements or a better understanding of the present data will 
clarify about the systematics. Nevertheless what is emerging from the above results 
is that the A He data are reaching a precision of the order of few percents. This fact 
requires a similar effort in the theoretical analysis, in order to reduce the uncertainty 
on the predictions at least at the same level of magnitude. In a previous paper |i~Tl 



wc 



performed a thoroughly analysis of all corrections to the proton/neutron conversion rates, 



) 



(a) u e + n — > e~ + p , (d) V e + p — > e + + n 

(b) e~ + p — > v e + n , (e) n — ► e" + 17 e + p , 

(c) e + + n^V e +p , (/) e~ +T7 e +p ^n , (1.7) 

which fix at the freeze out temperature ~ 1 MeV the neutron to proton density ratio. 
The Born rates, obtained in the tree level V — A limit and with infinite nucleon mass, 
have been corrected to take into account basically three classes of relevant effects: 

i) electromagnetic radiative corrections, which largely contribute to the rates of the 
fundamental processes, in particular in the low temperature regime, T < 0.1 MeV; 



ii) finite nucleoli mass corrections, which are of the order of T/M^ or m e /Mjv, with 
m e , M^v the electron and nucleon mass, respectively; 

iii) plasma effects, proportional to the surrounding plasma temperature, which both 
affect the microscopic process rates (a) — (/), as well as the neutrino to photon 
temperature ratio through e ± , 7 equations of state. 

The other main source of theoretical uncertainty comes from the partial knowledge of 
nuclear rates relevant for nuclei formation. Their numerical expressions, obtained by 
a convolution of the experimental data with a Boltzmann distribution, are affected by 
uncertainties of the order of 10% (see references quoted in [l^Q. More crucially, in many 



cases, these fits are known to well describe the data in a temperature interval which is only 
partially overlapping the one relevant for BBN, 0.01 MeV < T < 10 MeV. However, both 
a Montecarlo analysis to sample the error distribution of the reaction cross sections [[H| , 



and a more recent method based on linear error propagation [13], show that, in particular 
for 4 H e mass fraction, the effect is at most as large as the one due to the uncertainty on 
neutron lifetime r n , and smaller than 1%. Therefore it is theoretically justified to look, as 
in Jll| , for all sources of theoretical uncertainty up to this level of precision. The situation 
gets worse with D and Li, where the uncertainties due nuclear reactions can be as large 
as (10 -=-30)% 0. 



This paper represents the natural companion to fil l- We have built a new updated 
version of the standard BBN code, which is available since many years [[I], |P2| , where all 
corrections i)— iii) have been included. In particular we have also included the modified e ± , 
7 equations of state due to electromagnetic mass renormalization. In Section 2 we review 
the corrections to n ^ p Born rates, while in Section 4 we discuss the numerical method 
we have used to integrate the set of equations relevant for BBN, which are described in 
Section 3. The numerical results for light nuclei abundances, as functions of the final 
baryon to photon density ratio, 77, the number of effective neutrino degrees of freedom, 
N v , and the neutron lifetime, r n , are reported in Section 5, where they are discussed and 
compared with the experimental data. We have also performed a fit of these abundances 
with a precision of the order of 0.1% in the interesting range for the parameters rj, N u 



and T n . Finally in Section 6 we give our conclusions. 

2 Corrections to n <-> p Born rates 

As is well known, the key parameter in determining the primordial A He mass fraction, 
I4, is the value of the neutron to proton density ratio at the freeze-out temperature 
T ~ 1 MeV, since almost all residual neutrons are captured in 4 He nuclei due to its 
large binding energy per nucleon. In order to make an accurate theoretical prediction 
for Y4 it is necessary, though not sufficient, to have a reliable evaluation of the rates for 
the processes ( |1.7| ). An effort in this direction has been pursued in the last ten years by 
many authors. Recently, the entire set of corrections to the Born rates uob at the level 



of 1% accuracy have been recalculated in [] 1| and [15], with quite compatible results. In 
this Section we shortly summarize the main corrections Auj/ujb coming from considering 
radiative, finite nucleon mass, and thermal effects. This short review is here included for 
the sake of completeness and to fix the notation. A detailed discussion of the subject can 
be found in our paper [Tl[ . 



2.1 The Born rates 

Let us consider as an example, the thermal averaged rate per nucleon for the neutron 
decay process (e). In the simple V — A tree level, and in the limit of infinite nucleon mass, 
which we will refer to as Born approximation, one has 



uj B {n — > e +v e +p) 



& p ( C v 



3C A 



2tt 3 
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where Gp is the Fermi coupling constant, C v and C A the nucleon vector and axial coupling. 
In our notation p' and p' are the electron momentum and energy, and go the neutrino 
energy. The integration limits are imposed by the 0-function, go > 0. For reaction (e) 
we have go — M n — M p — p' Q = A — p' . The Fermi statistical distributions for e ± and 
neutrinos in the comoving frame, neglecting chemical potentials, are 



Fe( P ' ) 



JHl 



1 



Fv(q ) 



Algol 



+ 1 



(2.2) 



with (3 = 1/T and (3 V = l/T u []. All other rates for processes (a) — (d), (/) can be simply 
obtained from ( |2.1| ) properly changing the statistical factors and the expression for go [0 • 
The accuracy of Born approximation can be tested by comparing, for example, the 
prediction for neutron lifetime with the experimental value t^ x = (886.7=1=1.9) s [ |16|| . 
Using C v = 0.9751±0.0006 and C A /C V = 1.2601±0.0025 @, Eq. (Q in the vanishing 
density limit gives r n ~ 961 s. Therefore, to recover the experimental value, a correction 
of about 8% is expected to come from radiative and/or finite nucleon mass effects. In 
the same way these corrections are also expected to contribute to all six processes (o)- 
(/) relevant for BBN. In addition to these, microscopic n ^ p reactions taking place in 
the early universe, also feel the presence of the surrounding plasma of 7 and e^ pairs in 
thermodynamical equilibrium. Emission and absorption of real 7 or e from the thermal 
bath can be taken into account using the finite temperature field theory in the real time 



formalism. This has been considered by several authors |17|], and recently in [11 



2.2 Electromagnetic radiative corrections 

It is customary to separate the electromagnetic radiative corrections to the Born ampli- 
tudes for processes (|1.7|) in outer and inner terms. The first ones involve the nucleon 
as a whole and consist in a multiplicative factor to the modulus squared of transition 

amplitude of the form 

a 
1 + — <7frUo) • (2.3) 



The function g(p' ,qo) HH depends on electron and neutrino energies and describes the 



deformation in the electron spectrum. Its effect on a freely decaying neutron is such to 
reduce the Born prediction for the lifetime of about 1.6%. 

Inner corrections are sensible to nucleon structure details, and thus much more difficult 
to handle. They have been estimated in |TjJ, studying corrections to the quark weak 
currents. Translating the quark-based description in the hadronic language, the inner 
corrections result in the additional multiplicative factor 



-'^The ratio T v /T is fixed by entropy conservation and using the neutrino decoupling temperature flllj 



(see Section 2.4.1) 



where the first term is the short-distance contribution and A g = —0.34 is a perturbative 
QCD correction. The other two terms are related to the axial-induced contributions, with 
Ma = (400 -T- 1600) Me V a low energy cut-off in the short- distance part of the jW box 
diagram, and C related to the remaining long distance term. 

The global effect of these two kind of corrections, improved by resumming all leading 



logarithmic corrections a n ln n (Mz) [ 20| , is via the multiplicative factor 
G(Poi 9o) 



s( h £ + *0 + 2 £ i fc«-> + 'J 



S(M P ,M Z ) , (2.5) 



where a(/i) is the QED running coupling constant defined in the MS scheme and S(M P , Mz) 
a short distance rescaling factor, defined in [11] . 

Another effect to be considered, which can be in fact as large as few percents of the 
Born rates, is the so-called Coulomb correction, due to the rescattering of the electron in 
the field of the proton and leading to the Fermi function for Coulomb scattering 

HPo) - ( 1 + ot>K 1JL\ . (2.6) 

However, this effect is only present when both electron and proton are in either the initial 
or final states, namely it only corrects the amplitudes of processes (a), (b), (e) and (/). 

One may wonder if including the effects given by (|2.5|) and ( |2.6|) the theoretical pre- 
diction for neutron decay is now compatible with the experimental results. Evaluating 
numerically the integral over the phase space one finds r„ = 893.9 s, still at variance with 
the experiment. Even adding all known sub-leading effects the agreement does not really 
improve |H]. As in Ref. [IT] we take the point of view of rescaling all the rates ( |1 . 7| ) . 



after including finite nucleon mass corrections (see next section), by the constant factor 
1 + 5 T = t^/t^ = 1.008, which should be regarded as an energy independent correction 
to the weak process rates. This renormalization of the coupling guarantees the correct 
prediction for r n . 

In Fig. 1 we report the Born rates u>b for n <-> p processes, while in Fig. 2 we plot 
the corresponding radiative corrections, Aujr/ujb- Their effect is particularly large, up to 
~ 8%, at low temperature. 



2.3 Finite nucleon mass corrections 

There are three additional contributions to the n <-► p rates which appear when one relaxes 
the approximation of infinitely massive nucleons. The leading effects are proportional to 
m e /M^ or T/Mpj, which, in the temperature range relevant for BBN, can be as large as the 
radiative corrections considered in the previous Section. This has been first pointed out in 



22| and then also numerically evaluated in [11]. At order 1/Mjy there are new couplings 



appearing in the expression of the weak hadronic current, the larger one coming from the 
weak magnetic moments of nucleons 

j r = ^^ p (p)^(p-q7u n (q') , (2.7) 

where, from CVC, / 2 = V u d(fi p — /i n )/2 = l.81V u d- Both scalar and pseudoscalar contri- 
butions can be shown to be much smaller and negligible for the accuracy we are interested 
in. At the same order in inverse nucleon mass power it is also necessary to include the 
deformation of the allowed phase space for the relevant scattering and decay processes, 
due to nucleon recoiling. The sum of these two corrections for n <-> p rates with respect 
to th Born values, Aum/^b, is plotted in Fig. 3. 

The third effect is due to the initial nucleon thermal distribution. In the infinite 
nucleon mass limit, the average of weak rates over nucleon distribution is in fact trivial, 
since the nucleon is at rest in any frame. For finite M^, by considering only 1/Mjy 
terms, the effect of the thermal average over the thermal spreading of the nucleon velocity 
produces a purely kinetic correction Auk, whose expression can be reduced to a one 
dimensional integral over electron momentum which can be numerically evaluated. The 
explicit expression, which we do not report for brevity can be found in Section 4.2 and 



Appendix C of ]11|| . The ratios AujkI^b for n <-» p are reported in Fig. 4. Their size is 



rapidly growing with temperature, since they are proportional to the ratio T/Mn- 

2.4 Thermal-Radiative corrections 

The n <-> p rates, calculated as the processes would occur in vacuum, get slight corrections 
from the presence of the surrounding plasma of e^ pairs and 7. These are the so-called 
thermal-radiative effects. 

8 



To compute these corrections one may use the standard real time formalism for finite 
temperature field theory. According to this scheme, field propagators get additional con- 
tributions proportional to the number density of that particular specie in the surrounding 
medium. For 7 and e we have 

S(k 2 ) 



iS e (p') 



t^ + 2tt 



A- 2 
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o0\k o \ _ I 
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\xv 



— + 2ir8(k 2 )B(k 
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2n 5(p' 2 - m 2 e ) F e (p' ) (f + m e 



(2.8) 
(2.9) 



f — m e 
The entire set of thermal corrections Autr, at first order in its typical scale factor, i.e. 



aT/m e , have been computed by several authors [|17]] with quite different results. We have 



recently performed this lengthy calculation in |TT|, to which we refer for all details, and 
we have found a good agreement with the original result of first reference in [ 171 , namely 
that they contribute to correct the Born rates only for less than 1%. 

2.4.1 Radiative corrections on neutrino temperature 



By assuming a sharp neutrino decoupling at T& = 2.3 MeV [2~3|, the ratio T v /T can be 
evaluated using entropy conservation |T]J . This leads to the expression 



T 



1{xD)+21(xD) 
1 



1/3 



T <T D 
T>T D 



(2.10) 



with 



f2.11^ 



T{x) = f (y 2 + 2yxf 2 (Ay 2 + 8yx + 3a; 2 ) [exp(a; + y)±l] X dy . 

According to our notation x a = m^/T and x® = m^/Tr, with a = 7, e (+ or — in the 
above integrand is for fermions or bosons, respectively). Note that m^ and mf are the 
effective masses that photons and e^ acquire in the heat bath due to their interactions 
with the background plasma (see Appendix A for details). 

In [|11[] the neutrino temperature T v as a function of photon temperature T was eval- 

and 77ig — m e + aT 2 /m e . 



uated by using in (|2.10| ) the approximated expressions 



m 



R 



This simplified expression for T v has been used in all previous sections in order to obtain 
the Born rates and their corrections as a function of T only. The difference between the 
neutrino temperature evaluated with the correct renormalized masses ( |A.1| ), (|A.2| ), and 



the one obtained with the simplified expressions results to be smaller than 0.01%. The 
corresponding effect on the rates ( |1.7| ) due to this small change in neutrino/photon tem- 
perature ratio, which can be seen as a further sub-leading thermal-radiative correction to 
Born rates, can be neglected at the level of precision we are interested in. 

All thermal-radiative corrections to Born rates Aujtr/ojb are reported in Fig. 5 . 
As evident from this plot, around the freeze-out temperature T ~ 1 MeV, Auj^r only 
contribute as ~ 0.4% to the total rates. Thus they are clearly subdominant. Note that 
changing To in the range (2 -=- 3)MeV only affects T u /T for less than 0.2%. 

2.5 The total rates for n <-► p reactions 

In Fig. 6 we report the total corrections Auj to Born rates ujb- In order to use the 
corrected n ^ p rates uj = ujb + Auj in the BBN code, it is useful to fit their expressions 
as a function of the adimensional inverse photon temperature z = m e /T, 

1 13 

u{n^p) = — exp(-g„ p ,2) J^a, z~ l , 0.01 MeV < T < 10 MeV (2.12) 

T n 1=0 

J ^ exp (-q pn z) Y}li bi z- 1 0.1 MeV < T < 10 MeV 
u(p^n) = I r n 1 ,. ;/ „ /ni ,. T/ (2-13) 



0.01 MeV < T < 0.1 MeV 



with 



(2.14) 



a 2 = 0.456817-10 1 
a 5 = -0.583644-10 2 
a 8 = 0.371109-10 1 
an = 0.723091-10- 5 
q np = 0.340994 , 

6 3 = 0.109230-10 3 
h = -0.225830-10 1 
b 9 = -0.205210-10~ 3 
612 = -0.130731-10- 9 

(2.15) 

The fit has been obtained requiring that the fitting functions differ by less than 0.1% from 
the numerical values, while it is also a good approximation to consider a vanishing rate 
uj{p — > n) for T < 0.1 MeV, see Eq. ( |2.13| ), since it is a rapidly decreasing function with 
T^0. 
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a = 


: 1 


01 = 0.160615 


a 3 = 


: -0.401109-10 2 


a 4 = 0.137254-10 3 


a 6 = 


: 0.655938-10 2 


a 7 = -0.162185-10 2 


a 9 = 


: -0.378497 


do = 0.223840-10- 1 


O12 : 


= -0.462476- 10~ 4 


013 = 0.186287-10- 5 


b x = 


= 0.199695-10 2 


b 2 = -0.671993-10 2 


64 = 


= -0.29589M0 1 


65 = 0.40783M0 2 


b 7 = 


= 0.146751 


6 8 = -0.185408-10- : 


ho = 


= 0.158424-10- 5 


6 n = 0.369573-10- 6 


bi3 = 


= -0.329060-10- 9 


q pn = 2.89858 . 



3 The set of equations for BBN 

Let us consider N nuc species of nuclides, whose number densities, Xj = rii/n B , are nor- 
malized with respect to the total baryon density ub- The different nuclides are ordered 
in the following way: n, H, D, 3 H, 3 He, 4 He, 6 Li, 7 Li, 7 Be, ... (for the complete list 
see Ref. [12]). Denoting with R(t) the universe scale factor, the BBN set of equations as 
functions of R, tib, T, X iy and of the electron chemical potential e = p e /T readsQ 



R = \3M^ [Pt + Pe + Pu + Pb] 1/2 , (3.1) 



h B n R /247T 1/2 

r 2 , 



3 d = -\It72 [Pj + Pe + Pu + Pb} 1/2 , (3.2) 



n B R V M 



L (%,*.) = y^fM , (3-3) 

T = - I 3— [p 7 + p 7 + p e + p e + 6(T - T D ) (p„ + p^) + p B ] 

Xi = V iVi r H -« ' ,t\ - r«_ fc , ' /, = TiiXA . (3.5) 

fki I WiV*.! iVilJV,-! i J 

In the previous relations p and p denote the energy density and the pressure of an arbitrary 
particle specie. The function L(z,y) in (EO) is defined as 



1 roc / e v e -y \ 

L ( z ,y) = 7; dxxVx 2 -z 2 [ - - , (3.6) 

2 J z \ e x + eu e x + e y J 

i,j,k,l = l,..,N nuc , and the i-th nuclide, with charge and atomic number [Z^ A^), has 
mass Mi and mass excess (M u is the atomic mass unit) 

AM, = M t - AM U . (3.7) 

Moreover, in ( |3.5| ) we are considering in the sum a reaction between JVj nuclides of type z 

and Nj of type j which results in Ni nuclides of type I and N k of type k, with its reverse 
2 We are using natural units % = c = ks = 1. 
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reaction. The energy density and the pressure of baryons take the form 

3 



PB 
PB 



n B 



m u + J2 (am, + |r) Xj 



n B T J2 X j 



(3.8) 
(3.9) 



Eq. fl3.1| ) is easily recognized as the Friedmann equation where we have neglected for 
simplicity the cosmological constant. Eq. ( |3.2|) rules the scaling on ng, whereas (|3.3| ) 
states the neutrality of primordial plasma. From entropy conservation one gets ( p.4|) , and 
( |3.5|) are the Boltzmann equations for the N nuc nuclide number densities. Note that the 
presence of the B-function in fl3.4| ) is connected with neutrino decoupling at T = T B -f\ 

In the set of equations (|3.1| )- (|3~5| ) one can safely substitute Eq. (|3.3|) with an analogous 
relation, obtained expanding the l.h.s. of (|3~3|) with respect to <f) e , 

e x 



/•oo 

L(z,y)~y / dx x V ' x 1 — z 2 

J z 



(e* + l) 2 
In this case, Eq. ( |3.3|) provides an explicit expression for 



= yf- L (z) 



,(T,n B ,X. 



3)i 



a- m e \ n By -^ 



T J T 3 



(3.10) 



(3.11) 



The consistency of this approach has been tested by means of an iterative check. 

The set of equations (p.l[)-(p75|) can be transformed in a set of N nuc + 1 differential 
equations with the dimensionless variable z = m e /T as the evolution parameter. For 
numerical reasons, it is also better to turn the variable n B into the dimensionless quantity 
h = n B /T 3 , which varies more slowly with z than n B . In terms of these new variables the 
BBN set of equations becomes 



dh 

dz 



l-H{z,h,Xj) G(z,h,Xj 



3h 

z 



dXj 

dz 



G{z, h, Xj 



(3.12) 
(3.13) 



where the function G is 
G(z,h,Xj) = 



£ Q (4/3 a - z ?§f) + AQ(z D - z)p v + |/i £i Xj 



E a (p a + P«) + |e(« D - z) Pu + h Ej Xj\ H + hEj (zAMj + |) Tj 



(3.14) 



3 We have assumed that all neutrinos decouple at the same temperature Tq. Actually muon and tau 
neutrinos decouple at a slightly larger temperature of 3.5 MeV, but nevertheless our approximation is 
largely consistent with the required precision on A He yields. 
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In the previous equations z D = m e /T D , a = e, 7, and we have considered the dimensionless 
Hubble parameter H — H/m e , 



ill > Y ■■ i7im 



,, 10, Wj, 




3 M P z 2 



p 7 + p e + p u + hi zM u + Y, (zAMj + -) Xj 



1/2 



(3.15) 

and the quantities M u = M u /m e , AMj = AMj/m e , Tj = Tj/m e . Energy densities 
and pressures have been adimensionalized dividing by T 4 . In Eq.s fl3.12p and ( ^.13| ) we 
have neglected the terms containing the derivatives of chemical potential. In Appendix 
A we report the expressions for p a and p a evaluated taking into account the 7 and e ± 
electromagnetic mass renormalization. As already mentioned in the previous section this 
effect, changing the 7 and e^ equations of state, slightly modifies the T v jT ratio too. In 
order to speed up the numerical computation a fit of p Q and p a as functions of z has been 
performed and is also reported in Appendix A. 

The initial value for ( [3.12| ) is provided in terms of the final baryon to photon density 
ratio 77 according to the equation 



2C(3) 11 2C(3) 

7T 2 Vin 4 7T 2 



h m = ^mn= L i^ ■ (3.16) 

-t-Z A -77-Z 



The condition of Nuclear Statistical Equilibrium (NSE), which is satisfied with high ac- 
curacy at the initial temperature T in = 10 MeV, is then fixing the initial nuclide relative 
abundances. From NSE one gets, for an arbitrary i-th nuclide with gj internal degrees of 
freedom, 

A,(r ta ) = |(c(3) 1 /|)" ,_1 ^ (^) i,A - 1, ,,-.^X--exp{|} , (3.17) 
where Bi denotes the binding energy. 

4 Numerical Method 

The most critical part of the BBN code concerns the method of numerical resolution of 
the set of differential equations ( 3.12J) , (|3.13|) . In fact, since at high temperatures nuclear 



reactions proceed in both forward and reverse directions with almost equal rapidity, the 
r.h.s. of (|3.13| ) results to be a small difference of large numbers. When this occurs the 

13 



Table | 



1) 


n 


<— > 


P 


2) 


T 


<— > 


3 #e 


3) 


p + n 


<— > 


D + 7 


4) 


n + D 


<— > 


T + 7 


5) 


n + 3 He 


<— > 


4 #e + 7 


6) 


n + 6 Li 


<— > 


7 Li + 7 


7) 


n + 3 He 


<— > 


T + p 


8) 


n + 7 Be 


<— > 


7 Li + p 


9) 


n + 6 Li 


<— > 


4 #e + T 


10) 


n + 7 Be 


<— > 


4 #e + 4 #e 


11) 


V + D 


<— > 


3 #e + 7 


12) 


p + T 


<— > 


4 #e + 7 


13) 


p + 6 Li 


<— ► 


7 Be + 7 


14) 


p + 6 Li 


<— > 


4 Fe + 3 He 


15) 


p + 7 Li 


<— > 


A He + 4 #e 


16) 


D + A He 


<— > 


6 Li + 7 


17) 


T + A He 


<— > 


7 Li + 7 


18) 


3 He + A He 


<— > 


7 Be + 7 


19) 


D + D 


<— > 


3 ife + n 


20) 


D + D 


<— > 


T + p 


21) 


D + T 


<— > 


A He + n 


22) 


D + 3 He 


<— > 


A He + p 


23) 


3 He + 3 He 


<— > 


A He + p + p 


24) 


D + 7 Li 


■<— >■ 


A He + A He + n 


25) 


D + 7 Be 


<— > 


A He + A He + p 



Table 1: The reduced network of nuclear reactions. 



numerical problem is said stiff. As a consequence, the step size is limited more severely 
by the requirement of stability than by the accuracy of the numerical technique. In 
other words, to preserve integration stability it is required to use a shorter step size than 
what would be dictated by accuracy only. In order to manage the problem, we use a 
NAG routine implementing a method belonging to the class of Backward Differentiation 
Formulas (BDFs) |24| . This is quite a new approach for BBN codes. In fact the standard 
code |l], [12|] uses instead the implicit differentiating method (backward Euler scheme) |24| 
for writing the r.h.s. of ( |3.13|) and a Runge-Kutta solver. 

Few comments on the different numerical methods are in order. Let us consider the 
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differential equation 

dy(t) 



f(t,y(t)) • (4.1J 



dt 

In the Runge-Kutta methods the solution at ti+\ is completely determined by its value at 
ti (one-step methods), namely the solver has no memory. A different approach is provided 
by a wide class of numerical methods referred to as multistep methods like BDFs. Here, 
the values of the solution at t k (k = i,i — 1, ...,i — p), y(t k ) = y k , previously computed, 
and the unknown value y(ti+i) = yi+i, are interpolated by a polynomial, P(t; yi+i, y%, ■■■), 
in order to approximate the solution and its derivative. Substituting in the differential 

equation, 

dP 

-jr{U+i] Vi+i, yi,-)^f{U+i, Vi+i) , (4.2) 

one obtains a family of BDFs, 

(U+i-ti)f{t i+ i, s/i+i) ~ P(t%+i] Vi+i, Vi,--) ~ P(U, Vi+i, Vi,--) 

= a Q y i+1 + a 1 yi + ... . (4.3) 

Two methods can be used for solving the previous equation in the implicit case, a^ ^ 0: 
functional iteration and Newton's method. In the former case some initial guess is taken 
for y i+ i and refined by iteration. In the latter case, one linearizes Eq. ( f4.3| ) by expanding / 
around yi. The new point, yi+i, is then found by inverting a matrix, in a way similar to the 
backward Euler scheme. The NAG routine implements both methods and incorporates 
an error control test, which drives the step-size adjustment. 

The nuclear reaction network used in the code includes all the 88 reactions between 
the 26 nuclides present in the standard code [l], |l2j. We used the same nuclear rate 



data of the standard code, which are collected and updated in |^5j. In order to reduce 
the computation time one can also choose a reduced network, made of the 25 reactions 
between 9 nuclides listed in Table |l|. Using the complete network we have verified that 
the reduced one affects the abundances for no more than 0.01 %, while it greatly reduces 
the evaluation time. 
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5 Results on Light Element Abundances 

The reliable numerical code just discussed can now be used to study the effect of the 
different corrections ton<-^p Born rates on light elements abundances. By definition 

Y2 ~Y 2 ' Y3 ~x- 2 ' n - e, m] x ' Y7 ~x- 2 ■ (5 - 1} 

In the first two rows of Table ||] are shown the predictions for Y 2 , I3, I4 and Yj at 77 = 
5- 10 -10 , corresponding to the complete n <-> p rates, u>r t, and to the Born approximation, 
ujb 0- As it is clear from Table |2|, the main effect of the corrections, which results into 
the enhancement of the n — > p conversion rate, is to allow a smaller number of neutrons 
to survive till the onset of nucleosynthesis. This ends up in a smaller fraction of elements 
which fix neutrons with respect to hydrogen. 

The effects on light element yields due to the various corrections with respect to 
the Born predictions are also reported in Table |[ For all nuclides the pure radiative 
correction Aur provides the dominant contribution, while the finite nucleon mass effects, 
the kinetic and the thermal-radiative ones almost cancel each other. Finally the last 
row reports the further contribution due to the additional term required to recover the 



experimental neutron lifetime [11 



If we make use of the results of |14| to quantify the uncertainties coming from nuclear 
reaction processes, we observe that only for F 4 the radiative correction affects the Born 
result by an amount larger than the theoretical uncertainties, including nuclear reactions. 
For A He mass fraction in fact, the theoretical uncertainty due to nuclear reaction rates 
is estimated to be of the order of 0.1% and thus comparable with the uncertainty due to 
the experimental error on neutron lifetime. For D, 3 He and 7 Li the uncertainty due to 
the poor knowledge of nuclear reaction rates is estimated to be of the order of (10 -r- 30)% 



14]], thus completely covering any radiative/thermal correction onn<-»p rates. 

In Fig. 7 the predictions on Y4 are shown versus 77 for N v = 2,3,4 and for a 1 a 
variation of t^ x . The two experimental estimates for the primordial A He mass fraction, 
I4 and I4 , as horizontal bands, are also reported. Fig.s 8 and 9 show the predictions 



4 Note that, according to our notation, with ojb we denote the pure Born predictions fornHp rates 
without any constant rescaling of coupling to account for the experimental value of neutron lifetime (see 



Section |2J) 
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Table 





Y 2 


Y 3 


n 


Y 7 


^Tot 


0.3638-10' 4 


0.1175-10" 4 


0.2446 


0.2814-10" 9 


UJ B 


0.3727-KT 4 


0.1184-10- 4 


0.2550 


0.2873-10" 9 


Au r 


-2.3% 


-2.8% 


-3.8% 


-1.9% 


Au M 


.2% 


.1% 


.3% 


.2% 


Au>k 


.2% 


.1% 


.3% 


.2% 


Aujtr 


-.6% 


-.1% 


-.7% 


-.4% 


5 t ujt 


-.3% 


-.1% 


-.6% 


-.3% 



Table 2: The predictions on light element abundances obtained with the numerical code 
for i] = 5 • 10~ 10 and N v = 3. In the lower rows the effect of the various corrections is 
reported. 

for D and 7 Li abundances. Note that, due to the negligible variation of Y 2 and Y 7 on 
small T n changes, no splitting of predictions for 1 a variation of r^ x is present. 

A fit, up to 1% accuracy, of the relevant observables Y 2 , Y" 3 , Y4 and Y 7 as a function 
of x = log 10 (10 10 ?7), N v and r n has been performed. The following expressions have been 
obtained 



10 3 -F 2 



10 5 -F, 



^2 a i x% + a 5 (Nv — 3) exp j a,Q x + a 7 x 2 \ , (5.2) 

i=0 J 

4 

Y^ o,i x l + a 5 (N v - 3) exp{a 6 a;} , (5.3) 

8=0 

5 

10-F 4 = ^OiX 1 + a 6 (r - r ex ) + a 7 (N u - 3) + a 8 x(r - r ex ) + a 9 x(N u - 3) (5.4) 

10 9 -F 7 = Y,a i x i + a i {N u -3) + a 5 x{N u -3) exp j £> 5+i x* 1 , (5.5) 

.i=0 J u=i J 



where r P . 



.7 s and the values of the fit coefficients are reported in Table 131. 



Neutrino decoupling has been shown by many authors |26| to be a process which still 
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Table 



Coeff. 


io 3 -r 2 


io 5 -r 3 


io-y 4 


10 9 -F 7 


a 


0.4854 


3.325 


2.209 


0.5419 


ax 


0.2919 


0.1496 


0.5548 


-0.5981 


a-i 


-0.3516 


1.597 


-0.6491 


-1.914 


a 3 


0.5048 


-1.923 


0.7661 


4.521 


a 4 


-0.4269 


1.312 


-0.5366 


0.1587 


a 5 


0.7772-10- 1 


0.1782 


0.1614 


-0.3256 


a 6 


-4.397 


-1.705 


0.2059-10- 2 


-4.102 


a-j 


0.5925 


/ 


0.1300 


5.072 


a 8 


/ 


/ 


-0.4156-10" 4 


-1.209 


Og 


/ 


/ 


0.7433-10- 2 


-0.6269 



Table 3: Values of the fit coefficients (|5.2|) - (|5.5|) for light element abundances. 

takes place when e ± pairs annihilate. This implies that neutrinos are in fact slightly 
reheated during this annihilation process and their final distribution in momentum space 



shows an interesting non equilibrium shape. In Ref. |27[] it is estimated that the effect 
on Y4 due to the inclusion of this slight neutrino heating is very small, SY4 ~ 1.5 10~ 4 , 
in the whole range 10~ 10 < r] < 10~ 9 . We have included this constant correction to I4 
prediction. 

From the fit reported in Eq. ( |5.4|) it is easy to quantify the theoretical error on Y4. 
Since this is basically due to the uncertainty on r n we have 

AF 4 (a 6 + a s x)Ar n 



< 0.1% 



(5.6) 



y 4 ion 

In Fig. 10, our prediction for Y^ of Eq. ( |5.4j ) with N v = 3 and r n = 885.3 s, is compared 
with an analogous fit, Y^, performed in ||10|| . The agreement between the two expressions 
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obtained by independent codes is up to 1% in the relevant range for rj. 

In Fig. 11 we present the results of a likelihood analysis of the theoretical predictions 
obtained with the numerical code for the four combinations of experimental results: a) 
high D, low 4 He; b) high D, high 4 He; c) low D, low 4 He; d) low D, high 4 He, and 
using the low value for 7 Li abundance. In particular, we plot the product of gaussian 
distribution for D, 4 He and 7 Li centered around the measured values and with their 
corresponding experimental errors, 

L(N „ X) = exp{-<^-^> 2 }exp{-W<^-^> 2 } 

x expj-'^j-^} . ' (5.7) 

Notice that all functions have been normalized to unity in the maximum. As is clear 
from the plots, the analysis prefers the high value of D (plots a) and b)). In both cases 
the distributions are centered in the range x G 0.2 -j- 0.4, but at N v ~ 3 for low A He and 
N v ~ 3.5 for high A H e. For low D the compatibility with experimental data is worse. Note 
that c) and d) distributions have been multiplied by a factor of 25 and 100 respectively 
and centered in the range x G 0.6 -f- 0.8, and at N v ~ 2 for low 4 He and N u ~ 3 for 
high A H e. The better agreement at la of the data set a) and b) with the theoretical 
predictions is basically due to the effect of 7 Li data which corresponds to values for rj 
compatible with low D data of c) and d) at 2a only. It should be mentioned however 
that these results only take into account experimental errors, so that the confidence level 
regions in the N v — x plane would be broader by convoluting the considered distributions 
with the ones containing the theoretical error. 

6 Conclusions 

In this paper a detailed study of the effects on light element yields of the radiative, finite 
nucleon mass, thermal and plasma corrections to Born rates ( |1.7| ) has been carried out. 
The aim of such an analysis was to reduce the error on, basically, Y4 to less than 1%, which 
is motivated by the most recent experimental determinations for A He abundance. This 
accurate analysis has been performed using an update version of the BBN standard code 
fl], [12|]. A different numerical approach, based on BDF techniques has been implemented 
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to solve the stiff Boltzmann equations for nuclei densities. The numerical results for A H e 



mass fraction almost confirm the computation reported in Ref. [IU|, while the theoretical 
error, also including the propagation of uncertainties on nuclear processes, as estimated 
in 



14] , is of the order of 0.1%. Our analysis shows that the preferred experimental values 
are high value for D and low one for A H e, in which case the distribution is centered at 
x ~ 0.3 and N v ~ 3. 
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A Radiative corrections to e ± , 7 equations of state 

In an accurate description of the primordial plasma it is important to consider the elec- 
tromagnetic correction to the e^ and 7 equations of state induced by the e and 7 mass 
renormalizat ion . 

As well known, the photon renormalized mass, up to first order correction in the 



electromagnetic coupling constant a, reads [££ 

2 /(i 



m (z) ~ m e 



Z V 7T 



dx 



l + e a 



1/2 



(A.l) 



and for e [[L7j] 



mf(z,y) ~ m e < 1 + 







7T2; Z 



T + 



2^ 



.I'-' 



Z z + 



2vV 



log A 



dx 



l + e x 



(A.2) 



where z = m e /T, y = -E e /T and 



A(x,y,z) 



x 2 y 2 



— (z 2 + \fx r ^z 1 \/y 1 



(A.3) 



x 2 ?/ 2 — (z 2 — \Jx 2 — z 2 \Jy 2 — z 2 ) 

Note that in Eq.s ( |A.1| ) and (|A.2[) one can neglect the contribution of electron chemical 
potential, <fi e = /x e /T, due to its small value. 

By using ( |A.1|) and (|A.2| ) in the expressions of p 7 , p 7 , p e , p e one gets the latter 
quantities as functions of z only. Since the e^ and 7 energy densities and pressures have 
to be used in a BBN code, in order to speed up the computation one can fit these quantities 
as function of z and use these fits in the evolution equations. The fitted expressions for 
the dimensionless electron energy density and pressure, p e = Pe/T 4 and p e = p/T 4 , in the 
range z £ [0.05, 8.52] (p e = p e = for z > 8.52), result to be 

p e (z) = 1.145 + O^Sl-lO -1 * - 0.14543z 2 + 0.25507-10~V - 0.54168-10~V 



-3^5 



-5^6 



7^8 



- 0.11263-10^2 a - 0.29742-10-°z° + 0.38331-10~V + 0.45263-HT'z' 

+ 0.19241-10~V - 0.96597-10 _1 V° - 0.19505-10~ 10 z n - 0.14079-10" 12 ^ 12 , 



(A.4) 



p e (z) = (0.3786 + 0.19126-10 _1 z - 0.63895-10^ 2 + 0.32085-10~V - 0.48501-10" 
- 0.16611-10 3 ^ 5 + 0.82922-10'V + 0.79884-KTV - 0.60619-10" 6 / 



V 
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-6„9 nmnoi m-7„10 , n qck«/1 i n-8Ul\ .,-0.13145 z 2 



- 0.19568-10~°,2 y -0.1092M0~V u + 0.38564-10~V 1 J e 

(A.5) 

Moreover, in the considered temperature range, one can show that p 7 = p 7 /T 4 only 
varies between 0.6580 and 0.6573, and p 7 = p 7 /T 4 between 0.2193 and 0.2187. Thus, 
for simplicity, one can assume p 7 constant and equal to the average value 0.6577, and to 
0.2190 for p 7 . 
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Figure 1: The total Born rates, ujb, for n — > p (solid line) and p 
line). This notation is adopted hereafter. 
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Figure 2: The radiative corrections to Born rates, Aur/ub, for n ^ p transitions. 
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Figure 3: The finite nucleoli mass corrections to Born rates, Aujm/^b, for n 
tions. 
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Figure 4: The kinetic corrections to Born rates, Auk/ujb, for n <-► p transitions. 
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Figure 5: The thermal-radiative corrections to Born rates, Autr/ujb, for n 
tions. 
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Figure 6: The total corrections to Born rates for n <->■ p transitions. 
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Figure 7: The 4 He mass fraction, Y4, versus 77. The three solid lines are, from larger to 
lower values of Y4, the predictions corresponding to N v = 3 and t* x = 888.6 s, 886.7 s, 
884.8 s, respectively. Analogously, the dashed lines correspond to N v = 4 and the dotted 
ones to N v = 2. The dotted and dashed horizontal band are the experimental values of 
Ref.s [|| and ||, respectively, with la interval. 
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CM 




Figure 8: The quantity Y 2 versus rj is reported. The notation used is the same of Fig. 
8. Due to the negligible dependence of Y 2 on small variations of t^ x no splitting of lines 
is present. The horizontal bands dashed and dotted are the experimental values of Ref.s 

il- 
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Figure 9: The quantity Y? versus 77. The notation used is the same of Fig. 8. There is no 
splitting of lines related to Ar n , due to the negligible dependence of Y? on small variations 
of t* x - The horizontal bands dashed and dotted are the experimental values of Ref.s || [7| 

and ||, respectively. 
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Figure 10: The ratio (F 4 — Y[)/Y A versus log 10 (10 10 ^) for N v 
(see Section ET). 



3 and r„ = 885.3 s 10 
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Figure 11: The likelihood distributions for the light element yields Y2, Y4, Y7 are shown as 
functions of N v and log 10 (10 10 )?7, normalized to unity in correspondence of the experimen- 
tal values. From left to right and from top to bottom the following cases are considered: 
a) high D, low A He; b) high D, high 4 He; c) low D, low A He; d) low D, high 4 He. The 
plots for cases c) and d) are rescaled by a factor 25 and 100 times, respectively, compared 
to the one of a) and b). 
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